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Abstract. A direct numerical solution of the radiative transfer equation or any kinetic equation 
is typically expensive, since the radiative intensity depends on time, space and direction. An 
expansion in the direction variables yields an equivalent system of infinitely many moments. A 
fundamental problem is how to truncate the system. Various closures have been presented in 
the literature. We want to generally study moment closure within the framework of optimal 
prediction, a strategy to approximate the mean solution of a large system by a smaller system, 
for radiation moment systems. We apply this strategy to radiative transfer and show that 
several closures can be re-derived within this framework, such as Pjvi diffusion, and diffusion 
correction closures. In addition, the formalism gives rise to new parabolic systems, the reordered 
Pjv equations, that are similar to the simplified P]v equations. Furthermore, we propose a 
modification to existing closures. Although simple and with no extra cost, this newly derived 
crescendo diffusion yields better approximations in numerical tests. 



1. Introduction 



In many fields, macroscopic equations can be derived from mesoscopic kinetic equations. For 
instance, in the Navier-Stokes and Euler equations, the macroscopic fluid variables, e.g. density 
and momentum, are moments of the phase space distribution of the Boltzmann equation. Simi- 
larly, in the equations of radiative transfer |29j , the direction dependent kinetic equations can be 
transformed into a coupled system for the infinitely many moments. This can be interpreted as 
considering an (infinite) expansion in the kinetic variable and then taking only finitely many mem- 
bers of this expansion. This usually means that several "coordinates" of the field are neglected. 

A common feature of existing closure strategies is that they are based on truncating and approx- 
imating the moment equations, and observing to which extent the solution of the approximate 
system is close to the true solution. The approximations are supported by physical arguments, 
such as higher moments being small or adjusting instantaneously. In Sect. [2] we derive a moment 
model for radiative transfer, and in Sect. |3] we outline some classical linear closure strategies. 

A more potent strategy would be to have an identity for the evolution of the first N moments, and 
then deriving closures by approximating this identity. Moment models have first been studied from 
a higher perspective by Levermore [28] . Recent systematic studies of moment systems include also 
Struchtrup's order of magnitude method [35j which leads to the R13 equations of gas dynamics 
[36l[37]. In this paper, we take a similar approach. Extending results of [21], we show that the 
method of optimal prediction [TB] [THl HH [HI [S] can be applied to the equations of radiative transfer 
and yields closed systems of finitely many moments. Optimal prediction, outlined in Sect. [4] 
approximates the mean solution of a large system by a smaller system, by averaging the equations 
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with respect to an underlying probabihty measure. It can be understood as removing undesired 
modes, but in an averaged fashion, instead of merely neglecting them. 

Optimal prediction has been formulated for Hamiltonian partial differential equations [HI [T4] . 
however, without exploiting the full formalism. It has been applied to partial differential equa- 
tions |15l [2], however, only after reducing them to a system of ordinary differential equations 
using a Fourier expansion or a semi-discretization step. In addition, most considered examples 
are Hamiltonian, for which a canonical measure exists. In contrast, here we encounter partial 
differential equations (in space) after a Fourier expansion (in the angular variable). Hence, the 
methodology is generalized to semigroups. Furthermore, the radiation system is linear. Using 
Gaussian measures, the formalism is linear, and it yields an identity in the presence of a memory 
kernel. Here we restrict ourselves to Gaussian measures with vanishing covariance. We present 
linear optimal prediction in function spaces, and show various possible approximations of the 
memory term. 

In Sect. [5] we apply linear optimal prediction to the radiation moment system, and derive existing 
and propose new closure relations. The new formalism allows a better understanding of the errors 
done due to the truncation of the infinite system. The performance of the newly derived closures 
is investigated numerically in Sect. [6] 



2. Moment Models for Radiative Transfer 
The equations of radiative transfer [29] are 

-dti(x,n,t) + n\'i(x,n,t) + (a(x) + K(x))i(x,n,t) = f i(x,n',t)dn' + q(x,t) . (i) 

In these equations, the radiative intensity I(x, f2, t) can be viewed as the the number of particles at 
time t, position x, traveling into direction fi. Equation ([T]) is a mesoscopic phase space equation, 
modeling absorption and emission (K-term), scattering (cr-term), and containing a source term q. 
Due to the large number of unknowns, a direct numerical simulation of ([T]) is very costly. Often 
times only the lowest moments of the intensity with respect to the direction H. are of interest. 
Moment models attempt to approximate ([T]) by a coupled system of moments. 

For the sake of notational simplicity, we consider a slab geometry. However, all methods presented 
here can be easily generalized. Consider a plate that is finite along the x-axis and infinite in the 
y and z directions. The system is assumed to be invariant under translations in y and z and in 
rotations around the x-axis. In this case the radiative intensity can only depend on the scalar 
variable x and on the azimuthal angle 9 — arccos(/i) between the a;-axis and the direction of 
motion. Furthermore, we select units such that c = 1. The system becomes 

dtl{x,fi,t) + fidj{x,fi,t) + {(t{x) + K{x))I{x,fi,t) ^ J I{x,n',t)dfi' + q{x,t) (2) 

with t > 0, X G]a, b[, and /i e [—1, 1]. The system is supplied with boundary conditions that either 
prescribe ingoing characteristics or are periodic boundary conditions 

-^(a,M,^)U>o = f/(a,/i,i) = /(6,/i,i) 

-^(&,M,^)U<o = \lx{a,fi,t) = I^{b,fj.,t) 

and initial conditions 

l{x,fj.,0) = i{x,fj,). 

Under very general assumptions, this problem admits a unique solution / S C([0, t]; i^([a, &] x 
[— 1, 1])), i.e. at each time / is square-integrable in both spatial and angular variables [16j. 
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Macroscopic approximations can be derived using angular basis functions. Commonly used are 
Legendre polynomials |17l [7]. since they form an orthogonal basis on [—1,1]. Define the moments 



I{x,fj.,t)Pi{^i) d^. 



Multiplying ^ with Pk and integrating over fi gives 

dtlk{x,t) +dx J fJ'Pk{fJ')I{x,fJ.,t)d^i+ {a{x) + n{x))Ik{x,t) = {<j{x)Io{x,t) + 2q{x,t))5kfi. 
Using the recursion relation for the Legendre polynomials leads to 

oo 

dJk + dx^ bkili 



1=0 



k(2(?-/o), k^O 
-{K + a)Ik, k>0 



where bki — 2k+i ^k+i.i + iFkl'^*-'-!.'- This is an infinite tridiagonal system for k — 0,1,2, . 

dflk + bk^k-idxlk-i + bk.k+idxik+i = —Cklk + <Zfc , 
of first-order partial differential equations. Using the (infinite) matrix notation 

\ 



(3) 



B 



/O 1 

1 Q 2 

3 " 3 

i 

3 
7 



V 

we can write ([3| as 



C = 





(2Kq\ 







and q — 




J 





(4) 



+ B • = -c • I + q 



(5) 



The infinite moment system ([s]) is equivalent to the transfer equation Q , since we represented an 
function in terms of its Fourier components. In order to admit a numerical computation, (|3| 
has to be approximated by a system of finitely many moments Iq, . . . ,In, i-e. all modes /; with 
I > N are not considered. Mathematically, it is evident that such a truncation approximates the 
full system, since the /; decay faster than j, due to 

oo 

ll^(^,-,i)lli.(-M) = E^^'(^'^)'<°«- 



3. Moment Closure 



In order to obtain a closed system, in the equation for Ij^, the dependence on In+i has to be 
eliminated. A question of fundamental interest is how to close the moment system, i.e. by what 
to replace the dependence on In+i- In the following, we name three types of linear closure 
approaches: the P/v closure, higher-order diffusion (correction) approximations, and the simplified 
P/v closure. 



3.1. Pn closure. The simplest closure, the so-called Pn closure [6^ is to truncate the sequence 
i.e. /; = for / > N. The physical argument is that if the system is close to equilibrium, then the 
underlying particle distribution is uniquely determined by the lowest-order moments. This can be 
justified rigorously by an asymptotic analysis of Boltzmann's equation [5] . 
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3.2. Diffusion correction closures. The classical diffusion closure is defined for iV = 1. We 
assume Ii to be quasi-stationary and neglect /; for I > I, thus the equations read 

dtlo + dxh = -kIq + qo 
^d^Io = -(k + (t)/i . 

Solving the second equation for Ii and inserting it into the first equation yields the diffusion 
approximation 

dJo - dx^^^dxh = -nlo + go ■ (6) 

A new hierarchy of Pm approximations, denoted diffusion correction or modified diffusion closure, 
has recently been proposed by Levermore |28j . In slab geometry, it can be derived in the following 
way: We assume that /; = for ^ > + 1. Contrary to the P/v closure, the {N + l)-st moment 
is assumed to be quasi-stationary. Setting dtlN+i — yields the algebraic relation 

J — 1 N+l p, J 

which, substituted into the equation for In, yields an additional diffusion term for the last moment: 

T I -i a T /«(2'Z-/o), N = 

Ot^N ~T Om N — lOx^N—l — — I — V nOxx^ N — \ , , \' I 

XX n>Q ^ ' 

where — &jv,Af+i 2Af+3 ~ (2jv+i)'(2jv+3) ' ^^''^ ^ — ^ this closure becomes the classical diffusion 
closure (IgI. 



3.3. Simplified P/v closure. Other higher order diffusion approximations exist, such as the so- 
called simplified Pn (SPn) equations. For the steady case, these have been derived in an ad hoc 
fashion |2U [Ml lU and have subsequently been substantiated via asymptotic analysis [55] and via 
a variational approach |38l B|. In the time-dependent case, the simplified P/v equations have been 
derived by asymptotic analysis |20j . They are a system of parabolic partial differential equations. 
The SPi equations are identical to the diffusion approximation. The SP3 equations are 



3at 



. H- 202 - C] - <^a0 + q 



15a 

+ 202 



11 

210!' 



0-402 



12 



3a 



(l-a)-l C 



In these equations, a is a free parameter. To obtain a well-posed parabolic system, it should be 
set [20] to < a < 0.9. There are two particularly sensible choices: For a — ^, the equations 
reduce to the steady-state SP3 equations, whereas for a = |, they are asymptotically close to the 
time-dependent P3 equations. In the following, we set a to -k. 



These equations can be approximated by a quasi-steadiness assumption. This gives the simplified 
simplified P3 (SSP^) approximation 



dt(, 



Sat 



> + 202] - Cra0+ q 



2 

15a^ 



11 

21q?' 



3a 



O-t02 



3.4. Other types of closures. Various nonlinear approximations exist in the literature, most 
prominently flux-limited diffusion ^27J and minimum entropy methods |311 [TJ 1181 1391 119j . 
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4. Optimal Prediction 

Optimal prediction, as introduced by Chorin, Hald, Kast, Kupferman et al. [131 [lOj [121 |8j [9] seeks 
the mean solution of a time-dependent system, when only part of the initial data is known. It is 
assumed that a measure on the phase space is available. Fundamental interest lies in nonlinear 
systems, for which the mean solution decays to a thermodynamical equilibrium. Applications 
include molecular dynamics [34^ and finance |32] . The formalism has been developed in detail 
[TUl im [S] for dynamical systems 

x(i) = R(x(t)) , x(0) = X . (8) 

Let the vector of unknowns be split x = (x,x) into a part of interest x, for which the initial 
conditions x are known, and a part to be "averaged out" x, for which the initial conditions x are not 
known or not of relevance. In addition, let a probability measure /(x) be given (for Hamiltonian 
systems this could be the grand canonical distribution /(x) = exp{—(3H{x))). Given the 
known initial conditions x, the measure / induces a conditioned measure /j.(x) = Z~^/(x,x) for 
the remaining unknowns. An average of a function u(x, x) with respect to f~. is the conditional 
expectation 

, , f u(x, x) f (x, x) dx , , 

J /(x,x)dx 

It is an orthogonal projection with respect to the inner product (u,w) — E [uv], which is defined 
by the expectation E [u] = J J u(x, x)/(x, x) dx dx. Let (p{x,t) denote the solution of ([s]), for the 
initial conditions x = (x,x). Then optimal prediction seeks for the mean solution 

F(^(x,i) =E[^(x,i,t)|x] . (10) 



A possible, but computationally expensive approach to approximate ( 10 ) is by Monte-Carlo sam- 
pling, as presented in |12j . Optimal prediction formulates a smaller system for x. First order 
optimal prediction [TS' constructs this system by applying the conditional expectation ([9| to the 
original equation's (|8| right hand side. For Hamiltonian systems, the arising system is again 
Hamiltonian [9]. 

An approximate formula for the mean solution can be derived by applying the Mori-Zwanzig 
formalism [30' '40] in a version for conditional expectations [10 to the Liouville equation for 
([s]). It yields an integro-differential equation, that involves the first order right hand side, plus 
a memory kernel. First order optimal prediction can be interpreted as the crude approximation 
of dropping the memory term. Various better approximations have been presented in |111 1121 15]. 
For nonlinear systems, optimal prediction remains an approximation, even if the memory term is 
computed exactly. 



4.1. Linear Optimal Prediction. Consider a linear system of evolution equations (such as ([5])) 

dtu^Ru, u(0) = u, (11) 
where i? is a differential operator. Let the unknowns and the operator be split 



and R 



R R 
R R 



Lebesgue measures cannot be generalized to infinite dimensions. Instead, as in pSl [5]. we define 
measures on function spaces as follows. Let S* be a space of test functions (e.g. the Schwartz 
space) and let S' be the corresponding space of distributions, such that 5* C £^ C S" is a Gelfand 
triple jSj . Measures on vector valued function spaces are defined via a vector valued Gelfand triple 
5" C (L2)» c (S")". 
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Let f e S"" be a vector-valued test function, u e (S")" a vector-valued distribution, and S e 
(S")"^" a matrix-valued distribution. Then the pair-dot-product 



is a number, and the pair-matrix-vector-product 



is a vector. 



We define a Gaussian measure by its Fourier transform 

2 



y exp (i (f , u)^^,) dA(u) = exp ^— - 



+ i (f,m; 



SS' 



(12) 



where the variance S is a matrix-valued function and the expectation value m is a vector-valued 
distribution. 

Here we assume a measure with vanishing covariance between the different moments (i.e. diagonal 
S). As derived in [21^, this leads to the simple linear projections 

Px = Ex and (/ - P)x = Fx 

with 



E 



I 




and F 




I 



Considering the solution operators e*^ and e*-'^^, which we assume to be well posed, the Mori- 
Zwanzig formalism [SO] |40] yields the identity 

ft 



dt (e*«) = He 



tR 



tRF , 



RF+ / K{t~s)e'^ds 



(13) 
(14) 



where 

n = RE 

is the projected differential operator, and 

K{t) = e*-^^i?Fi?E (15) 

is a memory kernel for the dynamics. The projected operator e*^E represents the mean solution 
in the context of optimal prediction. Its evolution is given by 

ft 



dt (e*-"E) = TZe'^'E + f K{t ~ s)e"^Eds 

^0 



(16) 



Compared to (13 1, the middle term has canceled, since FE = 0. The projected differential 
operators, the projected right hand side (14 1 and the memory kernel (15 1 have the block form 

n 



HE, 



R 
R 



RE 



R 
R 



RERE 



RR 
RR 



and K 



K 
k 



Consequently, in ( 16 ) the evolution of the resolved moments is decoupled from the unresolved mo- 
ments. However, the resolved part of the memory kernel K{t) involves the orthogonal dynamics 
gtflF^ whose solution is as complex as the full system ( [TT] ). In the following, we derive new closures, 
by approximating the memory term in two steps. First, the integral in ( 16 1 is approximated using 
a quadrature rule. This is a general step that is independent of the specific form of the differen- 
tial operator R. We present various approximations in Sect. 4.2 Second, the exact orthogonal 
dynamics e*^^ are approximated by a simpler evolution operator. Good approximations employ 
the specific structure of RE. We present such an approximation in Sect. [5] 
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4.2. Approximations of the Memory Integral. The simplest approximation of the memory 
integral is to drop it. This yields the first order optimal prediction system 

dtu = 7^u , (17) 

which in some applications yields satisfying results. However, in most cases better results can 
be expected when approximating the memory term more accurately. This approach is generally 
denoted second order optimal prediction. Assume for the moment that the memory kernel K(t) 



be exactly accessible. The memory term in (16 1 uses the solution at all previous times. In 



principle, one could store the solution at all previous times, and thus approximate the integral 
very accurately (as presented in [33 J. However, such approaches are typically very costly in 
computation time and memory. In addition, for partial differential equations, the stability of 
such integro-differential equations (or delay equations, after approximating the integral), is a very 
delicate aspect. Therefore, here we focus on approximation that only use the current state of the 
system. We present three approximations. 

• Piecewise constant quadrature: Assume the memory kernel decays with a character- 
istic time scale t. Then a piecewise constant quadrature rule for K{s)u{t — s) yields 

K{s)u{t- s)ds K [ K{0)u{t)ds = TK{0)u{t) . (18) 
Jo 



• Truncated piecewise constant quadrature: For t < t, the memory integral in (18 1 
cannot reach over the full length r. Hence, it is reasonable to approximate 

t 

K{s)u{t - s) ds w min{r, t}K{0)u{t) . (19) 

• Modified piecewise constant quadrature: Using a piecewise constant quadrature rule 

for u{t — s) only, yields 

K{s)u{t'^ s)dsf=i / K{s)u{t)ds^ / K{s)dsu{t). (20) 
Jo Jo 



5. Application to the Radiation Moment System 



We now turn our attention to the infinite moment system (|3|. For consistency with the notation 



used in Sect. 



4.1 



we denote the (infinite) vector of moments by u(x, t) — {uq{x, t),ui{x,t), . . . 
In addition, we neglect the forcing term, since it is unaffected by any truncation. The radiation 
system (|3| can be written as 

dtU Ru , (21) 

where the differential operator R ~ —Bdx — C involves the (infinite) matrices We consider 
a splitting of the moments u = (u, u), where u — {uq, . . .) contains the moments resolved. The 
radiative intensity uq is always the starting moment. The other moments, however, could be 
reordered. The system sphts into blocks 



u = 


u 


, B = 


B 


B 


, and C = 


C 







u 




B 


B 





T 



1 



which we also assume to be the characteristic time scale of the system (r is a time 



where r 

scale since we have set c—1). The optimal prediction system (16 1 involves the block operators 



RE 



Bd^ + C 


, RF = - 




B5, 





Bdx 

Bdx + 



T J 



, and RFRE = 



BBd^x 
Md,, + ^Bdx 
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Here the upper left block of KE is the projected right hand side, which yields the first order 
optimal prediction system. The operator RF describes the orthogonal dynamics system 



dfU = B(?tU -u 



(22) 



Its solution operator e*^^^ is contained in the memory kernel (151. The equation for u in ( 22 1 is 
similar to the original full system (21), however, with one fundamental simplification. The decay 
matrix is a multiple of the identity, a fact that may be employed for fast approximations of the 
orthogonal dynamics. We shall pursue this approach in future work. For now, we construct a 
simple approximate system by neglecting the advection in the unresolved variables, i.e. B 0. 
The approximate orthogonal dynamics operator 



G 



RF 



has a simple solution 



u - r (1 - e-*/^) Bd^u 



Hence, the resolved component of the approximate memory kernel is given by 



e^'^RFKE = '^Bda;x - t (l - e"*/^) ^d^ (bBO^^ 
= e~*/^BBa^^ - T f 1 - e-*/^l BBBd^^a 



^Bd.. 



In the following, we work out the specific expressions for two different types of ordering of the 
moments. 



5.1. Standard ordering. We consider u = {uo,ui, . . . ,un) and u = (m^t+i, ■ ■ • )• With 

respect to this ordering, the projected right hand side 7?. is a simple truncated version of the full 

right hand side. The triple matrix product vanishes BBB = 0. And the double matrix product 
equals 

^0 ... 



BB 



... i9 



N J 



where i^n = (2N+'i){2N+3) • Hence, with the above simplified orthogonal dynamics, the optimal 



prediction integro-differential equation ( 16 1 equals 

ft 



dtu{t) = -Bda,u{t) - Cu{t) + [ e""/^BB9:,^u(i - s) ds . (23) 

^0 

Due to the special structure of BB, the memory term only modifies the A^-th equation, thus it is 
in fact a closure. The various approximations to the memory integral, presented in Sect. [4?2l yield 
the following modifications to the equation of the A'^-th moment: 

• First order optimal prediction: When dropping the integral term, the resulting sys- 
tem (23 1 becomes exactly the P^r closure. Note that by considering measures with non- 



vanishing covariances between the moments, all other linear closures can be derived 



Piecewise constant quadrature: The integral term in (23 1 is approximated by 



Ti9iv9.xwjv(t) . (24) 

A comparison with ([7| reveals that we have re-derived Levermore's diffusion correction 
closure [2H]. In the case A^ = 0, it is equivalent to the classical diffusion approximation. 
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Truncated/modified piecewise constant quadrature: The approximations (19 1, re- 
spectively ( pO| ) lead to the integral approximation 



fit)dNd^^UNit) . (25) 



Compared to the diffusion correction closure (24 1, the constant r is replaced by a time 
dependent function f(t), that increases gradually towards its maximum value r. The 
diffusion is turned on gradually over time. Hence, we denote these types of approximations 
crescendo diffusion {N — 0), respectively crescendo diffusion correction (N > 0) closure. 
More specifically, we call the approximation truncated crescendo diffusion (closure) when 



it is based on (19 1. In this case f(t) = min{r, t}. Similarly, the approximation based on 



(20), that leads to /(<) ~ t (l — e */^), is called modified crescendo diffusion (closure). 



The new crescendo closures ( |25[ ) introduce an explicit time dependence in the diffusion term. The 
physical rationale is that at t — 0, the state of the system (at least the unknowns of interest) 
is known exactly. Information is lost as time evolves, due to the approximation. This implies 
that crescendo diffusion closures shall only be used if the initial conditions are known with higher 
accuracy than the approximation error due to truncation. It is one way of circumventing an initial 
layer in time. In Sect. |6.1| we investigate the quality of these new approximations numerically. 



Although here we have assumed spatially homogeneous coefficients, we expect that the equations 
can be adapted to the space-dependent case in analogy to diffusion theory. Specifically, if ^(a;) 
and cr{x) are space dependent, we define t{x) = k{x)+it{x) ' replace tDxxUn by dx {T{x)dxUN)- 
The validity of this approximation will be addressed in future research. 



5.2. Reordered equations. In the following, we consider a reordering of the unknowns, given 
by the reordered vector v — Pu, where P is a permutation matrix that leaves — uq as the first 
component. This reordering is inspired by Gelbard's original derivation [2U[131[21] of the steady 
simplified Pn equations, in which the unknowns are the even-order moments and the odd-order 
moments are algebraically eliminated. The reordered system is 

dt'v ~ —Bdx'v — Cv , 

where we have the permuted advection matrix B — PBP"^. The decay term remains unchanged, 
since PCP^ — C. Of particular interest is the following even-odd reordering: For a given number 
N, consider the following permutation of the unknowns 

V = Pu = (Uo, "2, • ■ ■ , U2N, Ui,U3,..., U2N+l,U2N+2, • ■ ■ )^ > 

and consider the N + I lowest even moments v = (uq, U2, . . . , U2n) as of interest. The resulting 
reordered advection matrix has the following structure (here for N = 2): 









1 

2/5 3/5 
4/9 


5/9 








1/3 2/3 














3/7 


4/7 








B 






5/11 






6/11 










6/13 


7/13 














7/15 



The solid lines indicate separations between v and v. For the even-odd reordering, we have S = 0. 
As for standard ordering, the triple matrix product vanishes BBB — 0. The double matrix product 
D = BB becomes more interesting than for standard ordering. For any choice of iV, it is observed 
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to be positive definite. In addition, lower order matrices are contained as submatrices in higher- 
order ones. With the simplified orthogonal dynamics, the optimal prediction integro-difFerential 



equation (16) equals 



dtv{t) ^ -Cw{t) + f e-'*/^Da^^v(t - s) ds . (26) 
Jo 



The various approximations to the memory integral, presented in Sect. 4.2 now become: 



First order optimal prediction: Dropping the integral term in (26 1 yields a mere 
exponential decay for each moment 

dtv = -Cv , 

which is obviously a horrible approximation. For even-odd reordering, a non-trivial ap- 
proximation to the integral term is strictly required. 



• Piecewise constant quadratures: System (26 1 is approximated by 

atv = -6v + /(t)D5,,v, (27) 
where the coefficient function f(t) may be constant or time-dependent, depending on the 



integral approximation used. When based on (18 1, we have f{t) = r. In this case, we 



denote system (|27[) the reordered equations (RPn). When the approximation (19 1 



is used, we have f{t) — min{T, t}, and speak of the truncated crescendo RPn equations. 



Similarly, when (20 1, we have f{t) = r (l — e^*/"^), and speak of the modified crescendo 
RPn equations. In the case = 0, these systems become the corresponding diffusion 
closures (|6|. 



The RPn system ( |27[ ) has strong similarity with the time dependent SPn equations [20j. These 
are derived via asymptotic analysis. Although the unknowns are not directly comparable, a 
comparison of the respective diffusion matrices is compelling. For the choice of a = | the SP^i 
diffusion matrix is 

/0.3333 0.6667 -0.3333^ 
D = 0.1333 0.5238 
\0.3333 0.6667 0.2000 

In comparison, the RP2 matrix is 

/O.3333 0.6667 
D = 0.1333 0.5238 0.3429 
\ 0.1905 0.5065y 

Observe that the respective upper 2x2 submatrices are identical. The other entries differ, which 
can be expected, since the third unknown in the SP3 equations is not the fourth moment, as it 
is for the RP2 equations. Assuming the third variable in SP3 be quasi-steady leads to the SSP3 
system [20j, that turns out to be identical to the RPi system. Hence, similarly as both classical 
and new diffusion approximations can be derived from optimal prediction with standard ordering, 
even-odd ordering allows the derivation of existing and new parabolic systems that approximate 
the radiative transfer equations. In Sect. |6.2[ we investigate the quality of these new systems 
numerically. 



6. Numerical Results 



As a test case, we consider the ID slab geometry system ([2| on a: g]0,1[ (periodic boundary 
conditions) with k = a = 1.5, and no source q = 0. The initial conditions are uq{x,0) — 

exp (-500 ( X ~ k) ) and Uk{x,0) = Vfc > 1. The numerical solution is obtained using a second 

order staggered grid finite difference method with resolution Ax = 10^"^ and At = 0.8Ax. Any 
diffusive terms are implemented by two Crank-Nicolson half steps. 
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U(j(x) at t = 0.4 
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Figure 1. Evolution of energy distribution for various moment closures for = 



6.1. Diffusion Correction Approximations. The time evolution t e {0.1,0.2,0.3,0.4} of the 
radiative intensity uo{x, t) for various system sizes is shown in Fig. [ijfor = 0, Fig. [2]for iV = 1, 
and Fig.[5]for = 3. In each case, the thick solid curve is the solution to the full radiative transfer 
equation, obtained using a P/v closure, which has converged in N (TV = 51). The dotted function 
is the classical Pn closure (Sect. [3T| . The dashed curve is the diffusion correction closure ([t]), 
which for iV = equals the diffusion closure ([6]). The dash-dotted curve and the thin solid curve 



are the two versions of the crescendo diffusion correction (251, as derived in Sect. 5.1 The dash- 
dotted curve denotes the sharp time dependence, while the solid curve denotes the exponential 
time dependence. 

For iV = 0, the P/v closure gives a sole exponential decay of the initial conditions, which is a 
very bad approximation. The diffusion closure is significantly better than Pjv, but too diffusive. 
Both versions of crescendo diffusion remedy this problem fairly well, and yield good approxima- 
tions, especially for short times. Note that here the truncated crescendo diffusion reaches its full 
diffusivity at i = |. 

For > 1, the Pn solution splits into A^ peaks which move sideways with different velocities. This 
is due to the fact that the Pjv model is a linear system of hyperbolic equations with fixed wave 
speeds. With A^ increasing, the height of the peaks goes to zero, and the solutions approach the 
true solution. The diffusive closures smear out the highest considered moment u^, thus making it 
more uniform and reducing its influence on the second-highest moment u^-i- Hence, qualitatively 
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Figure 2. Evolution of energy distribution for various moment closures for iV = 1 
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Figure 3. Evolution of the 
error for moment closures for 
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Figure 4. Evolution of the 
error for moment closures for 
N=l 
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Up(x)att = 0.1 



true solution 
P., closure 



- diffusion corr. 
trunc. cresc. diff. corr. 
mod. cresc. diff. corr. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 5. Evolution of energy distribution for various moment closures for = 3 



an A^-th order diffusion correction solution lies "between" the Pjv-i and the Pjv solution. The 
results shown in Fig. [2] and Fig. [s] visualize this effect. In particular, one can observe that the 
classical diffusion correction is too close to Pn-i, i-e. the diffusion in the A^-th moment is too 
large. In contrast, the crescendo diffusion correction approximations ameliorate this effect. The 
peak height is reduced, and the approximate solutions are closer to the truth. 

The quality of the various approximations can be quantified by considering the time evolution of 
the L^-error of the approximate radiative intensity with respect to the true solution 

e(t) = r u;;pp""(:e, t) - t) dx . 

Jq 

The error of the various approximations is shown in Fig. [3] for A'^ = 0, in Fig. [4] for A^ = 1, and in 
Fig. [7] for A^ = 3. One can observe that the errors become smaller as A^ increases. Within each 
figure, the Pn closure yields the largest error at the final time t — 0.4. In contrast, the diffusion 
(correction) closures show smaller errors at t = 0.4. However, for short times, the error is increased 
compared to the Pn closure. This initial layer in time is particularly evident for A'^ = 0, shown in 
Fig.|3] 

Both versions of crescendo diffusion (correction) remove the initial layer of the diffusion closure, 
and the following increase in error is significantly smaller than with the Pn closure. In terms of 
accuracy, no significant difference can be observed between the two versions. 
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Figure 7. Evolution of the 
error for moment closures for 



Figure 8. Evolution of the 
error for parabolic approxima- 
tions 
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6.2. Parabolic Approximations. Fig.|6]shows a similar comparison for parabolic type approx- 
imations, as derived in Sect. |5.2[ The thick solid curve is the true solution. The dotted function 
is the classical diffusion approximation, which equals the RPq system. The dashed curve shows 



the RPi system (27 1, and the dash-dotted curve shows the RPg system. For increasing order, the 
reordered Pn equations converge to some limiting solution that is not equal to the true solution. 
While this is easily explained by the fact that the odd moments are not resolved, it implies that 
these systems are mainly useful for small values of N. For comparison, the thin solid curve shows 
the SP3 system (with a = ■^). One can observe that both the diffusion approximations and the 
SP3 equations are too diffusive. The RPn equations are far from perfect, but they do yield better 
approximations. 

The time evolution of the errors is shown in Fig. [s] For the considered examples, the errors 
obtained by the RPn equations are slightly smaller than the errors obtained by the diffusion 
closure or by the SP3. 



7. Conclusions and Outlook 



We have applied the method of optimal prediction to the infinite moment systems that model 
radiative transfer. An integro-differential identity for the evolution of a finite number of moments 
is obtained. The memory kernel involves an orthogonal dynamics evolution operator. Simple 
approximations of this operator and of the integral have been presented, and various approximate 
systems been derived. While traditionally closures had been derived using physical arguments 
or by asymptotic analysis, the optimal prediction formalism provides a very different strategy in 
which closures are derived from a mathematical identity. 

Using this methodology, existing closures can be re-derived, such as classical Pjy closures, diffusion 
and diffusion correction closures. In addition, new closures have been obtained. Fundamentally 
new aspects are crescendo diffusion correction closures, that modify classical diffusion correction 
closures by turning on the diffusion gradually with time. In the context of optimal prediction, 
this explicit time dependence has a natural interpretation as loss of information. While crescendo 
diffusion comes at no additional cost, numerical tests indicate that the quality of the results is 
generally improved. 

In addition, the application of the formalism to a reordered version of the moment system yields 
approximations of parabolic nature, which we denote reordered P/v equations. Some versions of 
the existing simplified P/v equations can be re-derived, and other systems can be obtained. A 
pleasing feature of of these new RPn equations is that the diffusion matrices are always positive 
definite. 

Future research directions will focus on better approximations of the orthogonal dynamics opera- 
tor, as well as of the memory integral. Approximations of the orthogonal dynamics can be based 
on the fact that higher order moments possess a uniform decay rate. The memory integral can be 
approximated using the solution at former times, thus leading to delay equations. While such sys- 
tems pose numerical challenges, it may be very rewarding to significantly improve approximations 
for very few considered moments. Another limitation of the presented analysis lies in its linear- 
ity. Nonlinear approximations, such as flux-limited diffusion |27j and minimum entropy closures 
[311 [TJ [TB] would require nonlinear projections. While the optimal prediction formalism is not 
restricted to linear projections, a precise formulation in the case of partial differential equations 
is yet to be presented. 
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